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ABSTRACT 

We calculate the chemical evolution of protoplanetary disks considering radial viscous accretion, 
vertical turbulent mixing and vertical disk winds. We study the effects on the disk chemical structure 
when different models for the formation of molecular hydrogen on dust grains are adopted. Our gas- 
phase chemistry is extracted from the UMIST Database for Astrochemistry (Rate06) to which we 
have added detailed gas-grain interactions. We use our chemical model results to generate synthetic 
near- and mid- infrared LTE line emission spectra and compare these with recent Spitzer observations. 

Our results show that if H2 formation on warm grains is taken into consideration, the H2O and 
OH abundances in the disk surface increase significantly. We find the radial accretion flow strongly 
influences the molecular abundances, with those in the cold midplane layers particularly affected. 
On the other hand, we show that diffusive turbulent mixing affects the disk chemistry in the warm 
molecular layers, influencing the line emission from the disk and subsequently improving agreement 
with observations. We find that NH3, CH3OH, C2H2 and sulphur-containing species are greatly 
enhanced by the inclusion of turbulent mixing. We demonstrate that disk winds potentially affect 
the disk chemistry and the resulting molecular line emission in a similar manner to that found when 
mixing is included. 

Subject headings: astrochemistry; protoplanetary disks; accretion, accretion disks; turbulence; in- 
frared: planetary systems 

1. INTRODUCTION 

Protoplanetary disks are a natural and active environ- 
ment for the creation of simple and complex molecules. 
Understanding their physical and chemical evolution is 
key to a deeper insight into the processes of star and 
planetary system formation and corresponding chemi- 
cal evolution. Protoplanetary disks are found encom- 
passing young stars in star-forming regions and typically 
dispe rse on timescales of less than a few million years 
re. g. JHaisch eTnil2?Ml20MlWatsoii et al.H2007lh Ob- 
servations of thermal dust emiss ion from these systems 
date back to the 198 0s (see, e. g.. lLada fc WillkingHl984l; 
iBeckwith et al.|[l990D and gave important in sight into the 
structure and evolution of disks (see, e. g.. iNattal I2007L 
for a review), however, inferring the properties of pro- 
toplanetary disks from dust continuum observations has 
several shortcomings: the dust component is only w 1% 
of the total mass contained in the disk and its proper- 
ties are expected to change with dust growth and planet 
formation. Also, dust spectral feat ures are too bro ad to 
study the disk kinematics (see, e. g.,[Carmona 2010( ). and 
observations suggest that the gas and dust temperatures 
can differ considerably, particularly in the UV - and X- 
ray heated surface layers of the inner disk (e. g.. lQi et all 
I20T31 ICarmona et al.ll200l . 
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The observation of molecular line emission from disks 
at (sub) mm wavelengths is limited by the low spatial res- 
olution and sensitivity of existing facilities and target the 
rotational transitions of abundant molecule s in the outer 
(> 5 AU) regions of the disk (see, e.g., iDutrev et al.1 
120071 for a review). Planet formatio n, however, mos t 
likely occurs at much smaller radii (|Armitagel 120071) . 
The superior sensitivity and spatial and spectral reso- 
lution of the Atacama Large Millimeter Array (ALMA) , 
currently under construction, is required to observe the 
(sub)mm emission lines from the inner regions of proto- 
planetary disks in nearby star-forming regions to enable 
direct study of the physical and chemical structure of the 
inner disk. 

Near- and mid-infrared observations, on the other 
hand, can probe the wa rm planet-forming regions (see, 
e.g., iNajita et all 120071 for a review). The 4.7/zm 
CO line emission has been detected at radii < 1 AU, 
using, for example, Kec k/NIRSPEC, Gemini/Phoenix 
and VLT/CRIRES (e. g ICardl2007t [Pontoppidan"et"aH 
[2008tlBrittain et alJ^OOa Ivan der Plas et al.ll2009t) . The 
Spitzer Space Telescope has been used to spectroscopi- 
cally study the inner disk at mid-infrared wavelengths. 
These observations show a rich spect rum of emission lines 
from simple organ i c molecules (e. g..[Carr fc Naiitall2008l ; 
iSalvk et al.1 120081: iPontoppidan et al.1 I2010D as well as 
from complex polycyclic aromatic hydrocar bons (PAHs) 
(e.g., lHabert et all [2001 iGeers et all l200l ) . They sug- 
gest high abundances of H2O, HCN, C2H2 and pro- 
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vide evidence that the disk supports an active chem- 
istry. C2H2 and HCN have also been detected in absorp- 
tion towards IR S 46 and GVTauri (jLahuis et al.l 120061: 
IGibb et al~ll2007l ). A particularly interesting case is CO2, 
which has been reported to be abundant in the circum- 
stella r disk of AATauri at ~ 1.2 AU (jCarr fe Najital 
120081) . in IRS 46 (lLahuis et al.ll2q06D an d in several other 
sources (e. g.. lPontoppidan et al.ll2010T ). but which is not 
detected in the majority of sources observed by Spitzer. 
The recently launched Herschel Space Observatory and 
future space telescopes (e. g., the Space Infrared telescope 
for Cosmology and Astrophysics (SPICA) and the James 
Webb Space Telescope (JWST)) will enable the observa- 
tion of molecular line emission in the mid- to far-infrared 
with unprecedented sensitivity and spectral resolution 
and will help address these controversial findings. 

Due to increasing computational power, the theoreti- 
cal study of the physical and chemical evolution of pro- 
toplanetary disks is entering its Golden Age. As a conse- 
quence of the observational limitations in the pre-Spitzcr 
era, many numerical models f ocused only o n the outer 
disk (e. g.. lAikawa et al.l [20021 : IWillacvl 120071 ) . With the 
launch of Spitzer, however, the che mistry of the inner 
disk h as increased in importance (e. g.. I Woods fe Willacy] 
I2009D . Unfortunately, complete models, including the 
physical and chemical evolution of the gas and the dust 
in the disk, coupled with a consistent treatment of the 
radiative transfer, is still beyond reach and several com- 
promises have to be made. 

One possibility is to ignore the chemical evolution 
and assume chemical equilibrium holds throughout the 
disk to investigate the physical structure in deta il 
(e.g., iKamp fc Dullemondl [20p: Uonkheid~el~aTl [20041) . 
Pioneering work bv iGorti fc Hollenbachl (|2004l ) and 
iNomura fc Millarl ([2005f i" combined a basic disk chemistry 
with a self-consistent description of the density and tem- 
perature structure taking into account the gas thermal 
balance, the stell ar irradiation a nd the dust grain prop- 
erties. Recently, iWoitke et al.l (l2009al ) calculated radi- 
ation thermo-chemical models of protoplanetary disks, 
including frequency-dependent 2D dust continuum ra- 
diative transfer, kinetic gas-phase and photochemistry, 
ice formation and non-LTE heating and cooling. These 
studies show, in accordance with the observations, that 
the stellar UV and X-ray irradiation increases the gas 
temperature considerably in the surface layer of the disk, 
altering the disk chemistry drastically in this region. 

Converse to that described above, several authors in- 
stead focused on combining the physical motion of mate- 
rial in the disk wi th time-dependent chemical network 
calcul ations (e. g.. lAikawa et al.l Il999t iMarkwick et al.l 
l2002Hllgner et al.ll2004l : INomura et alJl2009f ). Radial in- 
ward motion due to the viscous accretion and vertical 
motion due to inherent turbulent mixing or disk winds 
have the potential to change the disk chemistry signifi- 
cantly, depending on the timescales of the transport pro- 
cesses. Molecules, initially frozen out on dust grains, can 
evaporate and increase the gas-phase abundances, pro- 
vided that accretion and vertical mixing is comparable 
to, or faster than, the chemical reactions and the destruc- 
tion by incident stellar irradiation. First steps towards 
including the accre tion flow in the chem i cal network have 
been c arrie d out bvlAikawa et al.ldl999D.lM arkwic k et al.l 
(|2002T) and lllgner et al.l(|2004ft . INomura et al.l (|2009l ) cal- 



culatcd the effect of accretion along streamlines in steady 
1+1-dimcnsional a-disk models for chemical networks 
containing more than 200 species. They found that the 
chemical evolution is influenced globally by the inward 
accretion process, resulting in, e.g., enhanced met hanol 
abundances in the inner disk, m gner et all ()2004l ) also 
studied vertical mixing using the diffusion approximation 
and concluded that vertical mixing changes the abun- 
dances of the sulphur-bearing species. Their disk model, 
however, assumed the gas and dust temperatures were 
coupled everywh ere in the disk and ign ored stellar irradi- 
ation processes. iSemenov et al.1 (|2006[ ) investigated gas- 
phase CO abundances in a 2D disk model, allowing for 
radial and vertical mixing using a reduced chemical net- 
work of 250 species and 1150 chemical reactions, and find 
the mixing p rocesses strongly affe ct the CO abundances 
in the disk. iWillacv et "alT(|2006| ) studied the effects of 
vertical diffusive mixing driven by turbulence and found 
that this can greatly affect the column density of vari- 
ous molecules and increase the size of the intermediate, 
warm molecular layer in the disk. 

In this study, we consider the effects of inward accre- 
tion motion and vertical turbulent mixing to investigate 
to what extent the abundances of molecules such as H2O, 
CO2 and CH3OH are affected in the inner disk. As an 
alternative driving force for vertical transport, we inves- 
tigate the effect of upward disk winds on the disk chem- 
istry. Several kinds of disk wind m odels, such as mag- 
netocentrifugally driy en winds (e. g.. fBlandfor d fc Pavnd 
19821: iShu et al.ll 1994T ) and photoevaporation winds (e. g. 



Shu et aLlll993l) have bee n suggested thus far. Recently 



Suzuki fc Inutsukal (|2009l ) studied disk winds driven by 



MHD turbulence in the inner disk and found that large- 
scale channel flows develop most effectively at 1.5-2 disk 
scale heights, i.e., approximately in those layers where 
molecular line emiss ion is generated. The d isk model 
in our study is from INomura fc Millarl (120051) inclu ding 
X-ray heating as described in INomura et al.l (|2007f ) and 
accounts for decoupling of the gas and dust tempera- 
tures in the disk surface due to stellar irradiation. We 
use a large chemical network and allow for adsorption 
and desorption onto and from dust grains. We also cal- 
culate mid-IR dust continuum and molecular line emis- 
sion spectra assuming local thermodynamic equilibrium 
(LTE) and compare our results to recent observations. In 
Section [2] we describe our model in detail. In Section [3l 
we present and discuss our results, before drawing our 
final conclusions in Section [4] 

2. MODEL 
2.1. Disk structure 



Our physical disk model is from INomura fc Millarl 
(120051) with th e addi tion of X-ray heating as described in 
INomura et al.l (|2007l ). They compute the physical struc- 
ture of a steady, axisymmctric protoplanetary disk sur- 
rounding a typical TTauri star with mass M = 0.5M©, 
radius R = 2i? Q , effective temperature T e g = 4000 K 
and luminosity L = 0.87Lp . We adopt the a-disk model 
(jShakura fc Sunvaevlll973f l to obtain the surface density 
profile from an equation fo r angular momentum conser- 
vation (e.g.. iPringlel H98ll Eq. (3.9)). We assume an 
accretion rate of M = lO _8 M0yr -1 and a viscosity pa- 
rameter of a = 0.01. 
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We denote the radius in the cylindrical coordinate sys- 
tem with s to distinguish it from the radial distance from 
the central star, i.e., r — \/s 2 + z 2 , where z represents 
the height from the disk midplane. 

The disk is calculated between s = 0.01 AU and 300 AU 
in a 1 + 1-dimensional model, which corresponds to a to- 
tal disk mass of 2.4 • 1O -2 M for the values of a and M 
given above. The steady gas temperature and density 
distributions of the disk are obtained self-consistently by 
iteratively solving the equations for hydrostatic equilib- 
rium in the vertical direction and local thermal balance 
between heating and cooling of gas. To solve the eq uation 
for hydrostatic equilibrium (|Nomura et al.l [20071 Equa- 
tion (11)), the condition 



p gas (s,z)dz 



(1) 



with <p ga s(s,Zoo) 10 gem J is imposed 

(N omura fe Milled 120051) . The heating sources of 
the gas are grain photoelectric heating by FUV pho- 
tons and X-ray heating by hydrogen ionization, while 
cooling occurs through gas-gr ain collisions and by line 
transitions (for details, see INomura fe Millarl 120051 
Appendix A and INomura et al.ll2007T) . 

The dust temperature is computed assuming radia- 
tive equilibrium between the absorption of stellar radia- 
tion and resulting thermal emission by the dust grains. 
Viscous heating is only taken into account in the disk 
midplane since its effects in the disk surface are com- 
parably small for a = 0.01 (|Glassgold et all 12004) . In 
the midplane, the density is high and the gas and dust 
temperature are tightly coupled. Note that it is often 
assumed that the gas temperature is coupled with the 
dust temperature everywhere in the disk. This is not 
true in the disk surface which is strongly irradiated by 
the parent star. The dust temperature is calculated us- 
ing an iterative radiative transfer calculation by means 
of the short characteristic method in spherical coordi- 
nates, with initial dust temperature profiles obtained 
from the variable Eddington factor method (see iNomur, 
[200llDullemond fc Turollall2000l: iDullemond et alll200 
for details). 

The stellar X-ray and UV radiation fields are based on 
spectral observations of TW Hydrae, a classical T Tauri 
star. The X-ray spectrum is obtained by fitting a 
two-temperature thin thermal plasma model to archival 
XMM- Newton data (mekal model, L x ~ 10 30 ergs _1 for 
0.1 keV ... 10 kcV) and is slightly softer than that of typ- 
ical T Tauri stars with te mperatures of ab out 2 keV in 
their quiescent phase (c. f. JWolk et al.ll2005l ). The stellar 
UV radiation field has three components: photosphcric 
blackbody radiation, hydrogenic thermal bremsstrahlung 
radiation and strong Lya emission (Lfuv ~ 10 31 ergs _1 
for 6 eV ... 13 eV). The phot ochemistry u s ed in this pa- 
per is described in detail in iWalsh et all (|2010f ) . Note 
that the scattering of ly a emission by atomic hydrogen 
is not taken into account. On the other hand, the Lya 
emission does not affect the molecular abundances sig- 
nific antly when the du st grains are well mixed with the 
gas (|Fogel et aJ1l2011[j We inclu de the interstellar UV 
radiat ion field ( Draind (|1978l ) and Ivan Dishoeck fc Blackl 
(|1982l) l however, the interstellar UV and X-ray radia- 
tion fields are negligible compared to the central star 



(for details, see INomura et al.ll2007l) . 

The mass density, gas temperature and dust temper- 
ature for our disk model are displayed in Figure [1] top, 
middle and bottom, respectively. A logarithmic spac- 
ing of the grid cells is adopted in radial direction, while 
the vertical structure is calculated on a linear grid. The 
number of grid cells used in this work is 35 x 30 for 
0.5 AU < s < 300 AU and < z/s < 0.8. While this 
resolution is relatively coarse, its consequences for the 
line emission are limited due to the interpolation method 
during the line of sight integration (c. f., Section |2~7| . 

A more detailed description of the background theory 
and calculation of the disk p hysical model, alo n g wit h 
further references, is give n in |Nomura fc Millarl (|2005h . 
INomura et all (|2007T ) and IWalsh et all (|2010l) . ~ 



2.2. Dust properties 

The dust properties are important as they influence 
the physical and chemical structure of the disk in several 
ways. As the dominant source of opacity, they determine 
the dust temperature and the UV radiation field in the 
disk. They also affect the gas temperature, because pho- 
toelectric heating by UV photons is the dominant heating 
mechanism towards the disk surface. Further, the total 
surface area of the dust grains influences the molecu- 
lar hydrogen abundance through the H2 formation rates 
on grains and determines the freeze-out rate of chemical 
species onto the grain mantle throughout the disk (see 
Section [2.4p . Finally, due to their dominant contribution 
to the opacity, the dust grain model also determines the 
emerging continuum spectrum. 

For consistency, we adopt the same dust grain model 
for the spectral c alculations as for the u nderlying pro- 
toplanetary disk (jNomura fc Millarl [20051) . We assume 
the dust grains to be compact and spherical, and to con- 
sist of silicate grains, carbonaceous grains and water ice. 
We further assume that the dust and the gas are well 
mixed and that the gas-to-dust mass ratio is uniform 
(Pdust/Vgas = 0.01) in the disk. The size distribution is 
taken from iWeingartner fc Drainel (|2001a| ) to reproduce 
the observed extinction in dense clouds, with minimum 
and maximum sizes for the grains of 3.5 A and ~10 /xm, 
respectively. The heating of P AHs is included in the 
model for photoelectric heating (Weinga rtner fc Drainel 
2001bD. For further d etails, see also Appendix D in 
Nomura fc Millarl pOOBT l. 



2.3. Gas-phase reaction network 

For the calculation of the chemical evolution in the 
disk, we use a subset o f the UMIST Databa se for Astro- 
chemistry, or Rate06 (jWoodall et all l2007| ) , neglecting 
only those species, an d thus reactions, containing flu- 
orine and phosphorus (W alsh et al.ll2010l ). Since we are 
primarily interested in those regions where molecular line 
emission originates and where most molecules are in the 
gas phase, we neglect grain-surface reactions excepting 
H2 grain formation (see Section l2.4p . In summary, our 
gas-phase reaction network contains 375 atomic, molecu- 
lar and ionic species, ranging from atomic hydrogen H to 
relatively complex molecules such as H3CgN + with mass 
m = 125 u (see Table Q] for a list of species of particular 
interest in this study and their binding energies on dust 
grains). The reaction network consists of 4346 reactions 
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Fig. 1. — From top to bottom: Gas density, gas temperature and dust temperature for the full disk model (left) and for constant 
streamlines z/h of the viscous accretion (right), where h denotes the scale height of the disk at radius s. For the gas temperature, inlaid 
are the vertical profiles as function of z/s in the inner disk region. 



including 3957 two-body reactions, 214 photochemical 
reactions, 164 X-ray/cosmic-ray induced photoreactions 
and 11 reactions for direct cosmic-ray ionization. 

2.4. Gas-grain interactions 

The formation of molecular hydrogen on the surface 
of dust grains is thought to be the most efficient means 
of forming H2 directly from atomic hydrogen. Except 
for this process, we neglect grain-surface reactions of icy 



mantle species, but allow for adsorption onto and des- 
orption from dust grains. Hence, except for atomic hy- 
drogen, a particle that hits a dust grain can either detach 
without change or freeze out and remain in the icy grain 
mantle from where it may eventually evaporate. Atomic 
hydrogen forms H2, freezes out onto the grain mantle or 
evaporates without reaction. 

For the H2 grain formation rates, we consider two 
different models. In model A, we refer to early work 
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bv iGould k Salpeterl (|1963f) and iHollenbach k Salpeterl 
(|1971f ). who derive a grain formation rate for imperfectly 
shaped grain surfaces of 

n(H 2 )| r = in H ■ «H h ( T ga S )«d us t7r(a 2 )Q;H 



n H -2-10 



48.-1 / -*f 

S Tltot 



V300K 



1/2 



(2) 



In Equation {2}, ^H^gas) i s the thermal relative velocity 
of hydrogen atoms, 



( 8 



7rmH 



1/2 



(3) 



and /idust7r(a 2 ) the cross section of the dust grains per 
unit volume, with n du st ~ 6 ■ 10~ 10 nt o t and (a 2 ) = 
4.84 • 10 _12 cm 2 for the dust grain model adopted here. 
Further, n to t and «h denote the total number density of 
hydrogen nuclei and the number density of atomic hydro- 
gen, respectively. The total number density of hydrogen 
nuclei is calculated from the mass density, p gas , using 

Pgas Pgas 



ntot = 



(1 + 4y H c)"i H I Ami 



where j/ho = 0.1 is the fractional elemental abundance 
of helium. Here, we neglect the mass contained in other, 
less abundant elements. 

Alternatively, in mo del B, the H2 grain form ation is 
calculated according to iCazaux k Tielensl (|2002l ) , 



"(H 2 )| r : 



= n(H 2 )|* 



e(Tdust) 



= ^"h • v t ^{T g&s )n AvLSt Ti{a 2 )aYi£(T AvLSt ) , (4) 

with several updates for the efficiency e(Td us t) (Woitkc 
2009, private com munication). The sticking c oefficient, 
«h, is taken from IHollenbach k McKeel (|1979[ ). and ap- 
proximated to first order, i. e., 



1 + o.V(r ga8 + r dUBt )/iooK) 



(5) 



The model for adsorption of gaseous species onto dust 
grains and de sorption of ic e speci es from dust grains 
is taken from iWoitke et al.l ()2009aD , who consider col- 
lisional adsorption together with thermal desorption, 
cosmic-ray induced desorption and UV photodesorption. 

We calculate the ice formation rates for all neutral and 
negatively charged species k in our reaction network, 



„ pads 

= n k K k 



l k ice 



(6) 



R 



dcs, th 



R 



dcs, cr 



R 



dcs, ph 



where nfci CG denotes the number density of ice species k 
and 



-.desorb 



"k ice ^s fraction located in the uppermost surface 
5 of tl 
(jAikawa et al.lll996l) 



layers of the ice man tle. This fraction is given by 



„ desorb 



W'k ice 

W'k ice 



"act" 



Mcc ^ ^act 
Mce ^ ^act : 



(7) 



with n ac t = 3 ■ 10 15 • 47r(a 2 )?idust being the number den- 
sity of active spots in the ice mantle, and n- lce being the 



TABLE 1 

Molecular binding energies E bmd [K] 



species 



J [K] 



reference 



H 




350 


(5) 


H 2 




450 


(5) 


CO2 




2690 




H 2 




4820 


:i: 


CO 




960 


(i) 


OH 




1260 




CH 4 




1080 


8! 


NH 3 




3080 




HON 




4170 




CH3OH 




2140 




C 2 H 2 




2400 


(2) 


SO2 




3070 


(5) 


H 2 S 




1800 


(5) 


References. 


— (1) Sandford & Allamondola 


(1988), (2) 


Yamamoto et al. (1983), (3) Sandford & Allamondola (1990), (4) 
ISandford k. Allamondolal 119931). (5)IHaseeawa & Hcrbst Jl993) 



number density of all ice species (n to t = «gas + ^ice)- ln 
the following, we summarize the e xpressions for t he rate 
coefficients and we refer readers to IWoitke et al.l (|2009aD 
for a detailed discussion. 

The adsorption rate for species k with mass m k and 
thermal velocity wt h is 



Rf s = it(a 2 )n dust a k vl 



th 



(8) 



where we assume a sticking coefficient of unity for all 
species excepting hydrogen (c. f., an in Equation (|5j|). 
The desorption rates are 



R 



dcs, th 



:^ sc exp 



^bind 



T, 



dust 



R, 



hoxR 



fc,70K 



5 ■ 10" 



R 



dcs. ph / 2\ n dust 

= 7r(a ) 

"act 



^fc-pFUV • 



(9) 
(10) 
(11) 



The binding energies El ind [K] , for the most important 
species are listed in Table [U The vibrational frequency 
of species k in the surface potential well is 



3 • 10 15 fc B £]? indx 1/2 



TT 2 m k 



(12) 



The expressio n for cosmic-ray in d uced desorption is 
adopted from Hascga wa k Herbstl (|1993l ). who evalu- 
ate the rate coefficient at a temperature of 70 K with 
a duty cycle of the grains of /70 = 3.16 • 10 -19 . The 
photo-desorp tion yields Y k are s imply set to 10 -3 per 
UV photon (jWestlev et all 119951) . Finally, the cosmic 
ray ionization rate £cr and the FUV field -Ffuv in 
units [photons/cm 2 /s] are calculated from the density 
profile and the dust opacity adopted in the disk model 
(INomura et"aT1l2007D. The calcula tion of the FUV field is 
based on lNomura fc Milarl (|2005l . Equations (6) and (7)) 
and accounts for dust scattering. 

2.5. Initial conditions 

The initial abundances we use as input are listed in Ta- 
blc[2]and we calculate the chemical evolution over 10 6 yr, 
the typical lifetime of protoplanetary disks. The initial 
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TABLE 2 
Initial fractional abundances in 
the chemical network, including 
the gas and ice component of each 

SPECIES. 



species 


a 

x k,0 


H 


2.00c-07 


H+ 


1.00e-10 


H+ 




He+ 


2.50e-ll 


He 


1.00e-01 


CH 4 


2.00e-07 


NH 3 


1.00e-05 


H 2 


1.00e-05 


C2H2 


5.00e-07 


Si 


3.60e-08 


N 2 


2.00e-05 


C 2 H 4 


1.00e-08 


CO 


5.00e-05 


H 2 CO 


4.00e-08 


CH3OH 


5.00e-07 


2 


1.00e-06 


H 2 S 


1.00e-06 


Fe+ 


2.40e-08 



a The initial fractional abundance is de- 
fined as Xfc o = n k o/^tot- 



include both the ice 



abundances for each species rifc^ice. 
and the gas component. 

We compare the freeze-out timescale r^ ds and evapo- 
ration timescale t^ cs for each species in a simple model 
to determine whether it exists initially as gas or ice only: 



ads 



= R 



;u Is 



_dcs 



= R 



des^ " 



(13) 



For r^ ds < r d 



the species is initially frozen out, n k i cc — 
rife and n k = 0. Otherwise, it exists in the gas phase 
only, n k = n kfi and n kicc = 0. 

2.6. Physicochemical interactions: viscous accretion, 
turbulent mixing and disk winds 

Modeling the physicochemical interaction is the key 
aspect of our study. Even with the steady increase in 
computational power, this is still a challenging problem 
which can only be overcome by massive parallel comput- 
ing and by accepting several simplifications. In the con- 
text of a steady disk model, for example, it is a reasonable 
assumption to describe the interaction between the disk 
physics and chemistry as a one-way process. Hence, we 
consider the influence of physical processes on the disk 
chemistry, but neglect the feedback of the chemistry on 
the physical disk structure. 

Despite this simplification, the calculation of the disk's 
chemical evolution with interaction is still considerably 
more computationally expensive, the reason being that 
adjacent grid cells are chemically coupled through the ra- 
dial and vertical exchange of material. Accordingly, their 
evolution must be calculated in parallel with integration 
time steps that can be shorter than the the time scale 
of the chemical reaction network, due to the Cqurant - 
Friedrichs-Lewy (CFL) condit ion (jCourant et al.lll928| ) 
of the fluid (see Section I2.6.4p . More precisely, to com- 
pute the models with interaction, we adopt the following 
procedure: all grid cells are initialized at the starting 
time with their initial conditions (Section I2.5[) and the 
transfer coefficients between the grid cells are evaluated 



(see below for the calculation of these coefficients for the 
different transport processes). The chemical evolution of 
the entire disk is calculated in parallel for time steps that 
satisfy the CFL condition and the transport coefficients 
are updated each iteration. This step is repeated until a 
total age of the disk of 10 6 yr is reached. 

In the following, we describe each process we consider 
in this study. A detailed discussion of the numerical 
methods used to compute the chemical evolution is given 
in Appendix [AJ 

2.6.1. Radial transport by viscous accretion 

Firstly, we allow for the radial inward transport of 
material through the disk via accretion. The main pa- 
rameter describing this effect is the accretion rate M = 
lO _8 M0yr _1 , which relates the radial inflow velocity v s 
with the column density £ at radius s, 



M = 2irsv s Il (v s > for radial inflow) 



(14) 



We assume that inward accretion occurs along stream- 
lines I = Vs 2 + z 2 , with z/h = const. Here, h is the scale 
height of the disk defined by h = c s /Qk, c s is the sound 
speed at the disk midplane and fix is the Keplerian an- 
gular frequency 

We denote the net accretion rate of species k, i.e., the 
net change in its number density, with n kacc . In steady 
state, the net accretion rate of species k can be written 
as (see Appendix I A. 1 1 for the numerical implementation) 



rik, 



d{An k vi) 
di 



(15) 



where A is a geometrical factor and corresponds to 
Pi+i/Pi hi (jA2[) . The boundary conditions for the ac- 
cretion process at the inner and outer disk radii are set 
as follows: due to the purely inward motion, all ma- 
terial that flows inwards from the innermost grid cell 
is removed from the disk, whereas at the outer radius 
s = 300 AU, the accretion timescales become sufficiently 
long that we can safely assume n k _ acc = for all species. 
In Section 12.6.41 we compare the timescales for viscous 
accretion with the chemical timescale and also with cor- 
responding values for the other physical processes con- 
sidered here. 

2.6.2. Vertical transport by turbulent mixing 

Although the mechanism is not yet fully clear, it is 
commonly understood that the source of viscosity in as- 
trophysical disks is turbulence. Hence, a sufficiently large 
viscosity, as is the case in our model (a = 0.01) implies 
that the effect of turbulent mixing is strong, particularly 
in the inner disk region. While this turbulent mixing, 
in principle, occurs in the radial and vertical directions, 
we ignore the radial turbulent motion and focus on the 
vertical mass fluxes only. Furthermore, we assume that 
the turbulent velocity is independent of the distance z 
from the midplane and thus only a function of radius. 

We calculate the mixing in a diffusion-type approach. 
The diffusion coefficient is set as K = v z dz, where 

V z = Wz.mix = Ct'c s , (16) 

in analogy to the expression for the turbulent velocity 
in a Shakura-Sunyacv disk model, where the turbulent 
viscosity is 
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^turb — ^turb'turb 



factor x = 10 in Equation (Tl9[) to reduce the wind ve 
(17) locities such that the dispersal time scales, £/(/°g 



We assume a' = 10 ~ 3 to account for the turbulence in 
vertical direction only. For our numerical model, the 
length scale dz corresponds to the vertical extent of the 
grid cells Az as h/30. The mixing timescales are dis- 
cussed in Section [2.6.41 

The usu al diffusion equation for species k (c. f., Equa- 
tion (3) in lWillacv et~atll2006| ) reduces to 



d ( ts d{n k /n tot ) 

«.fc, mix = TT~ Kn tot 



(18) 



In Appendix IA.21 we describe the numerical implemen- 
tation of calculating the chemistry accounting for the 
mixing rates. 

The boundary conditions are set at the disk midplane 
and surface. Due to the symmetry with respect to the 
midplane, the net exchange rate at z — is zero. Like- 
wise, we assume that there is no exchange of material 
between the uppermost disk layer and any potential at- 
mosphere on top of it . 

2.6.3. Vertical transport by disk winds 

We consider vertical disk winds as an alternative de- 
scription of the vertical motion in the chemical model 
of the disk. Contrary to turbulent mixing, which im- 
plies upward and downward motion, the disk wind trans- 
ports mater ial upwards only (c. f. Figure [TlTc) in Ap- 
pendix E]2|). 

We refer to the wind solutions of iSuzuki fc Inutsukal 
(|2009f ) to describe the vertical disk wind in our protoplan- 
etary disk model. They investigated disk winds driven by 
MHD turbulence in local three-dimensional MHD simu- 
lations of stratified accretion disks and found that large- 
scale channel flows develop most effectively at 1.5-2 disk 
scale heights, where the magnetic pressure is comparable 
to the gas pressure (here, we refer to their model with 
initial plasma /3 = 10 6 ). The breakup of these channel 
flows drives structured disk winds, which might play an 
essential role in the dynamical dispers a l of p rotoplan- 
etary disks. While ISuzuki &: Inutsukal (|2009f l perform 
local three-dimensional simulations at radius s = 1 AU 
which lead to quasi-periodic cycles of 5-10 rotations and 
strong mass outflows, our disk model is steady and does 
not allow for a dispersal of the disk. 

We therefore adopt their solutions for the wind veloc- 
ity {>wmd,o at the midplane at 1 AU a nd calculate the ra- 
dial a nd vertical profile according to ISuzuki fc Inutsukal 
(|2009f ) and the vertical continuity equation in steady 
state, d(pv z )/dz = 0, 



Vo = V wilL d{s,Z = 0) = KV windt0 

p(s,z = 



ld (s,z) = Vq 



' 0,(1 AU)' 
0) 



p(s,z) 



(19) 



wher e Wwind.o = 

8.5 • 10 _5 c s (l AU) (|Suzuki fc Inutsukal 
l2009f h For large values of v w ind,0: the midplane layers 
disperse over timescales shorter than the typical disk 
lifetime of 10 6 yr, which contradicts the assumption of 
a steady disk model. We therefore introduce an ad-hoc 



(where £ is the surface density of the disk), become 
~ 10 6 y r. The resulting wi n d pro files reproduce the re- 
sults of ISuzuki fc Inutsukal (|2009l ) in the sense that the 
disk wind is most effective at ~ 2 disk scale heights, 
where the gas density and therefore the gas pressure 
decreases steeply and any potential magnetic pressure 
would become comparable with it . 

The wind rates are calculated from the vertical con- 
tinuity equation for species k, using the wind velocity 
profile (Equation (fT9jO : 



rik, 



d{n k v z ) 
dz 



(20) 



The wind rates for the discrete grid of the disk are given 
in Appendix IA.3I 

The boundary conditions are set as follows: since we 
adopt small wind velocities so that the disk dispersal 
time becomes long enough, we can assume that the mass 
loss is negligible at the midplane and therefore h k _ w ; n( j = 
0, whereas at the upper boundary, we simply allow all 
material transported by the disk wind to be lost. 

2.6.4. Timescales 

The grid crossing times for the mass transport by vis- 
cous accretion, turbulent mixing and disk winds are dis- 
played in Figure [H These are set by the CFL-condition 
for the discrete grid of our disk model: 



As 



Az 



Twind 



"Z. wind 



The radial and vertical size of the grid cells are given by 
As and Az, respectively, and depend on distance s from 
the star. For reference, we display the accretion, mixing 
and wind velocities in Figure [2] (left panel). 

In the right panel of Figure JZl we also plot the timesc ale 
of photochemistry given by rtPulev fc Williamdll98l : 

Tphotochem = 10 yr • — , (21) 

F FU v(s,z) 

where -Ffuv.isrf is the FUV photon flux of the interstel- 
lar radiation field. In the upper layers of the disk, the UV 
radiation field is strong and the timescale for the photo- 
chemistry is shorter than the dynamical timescales. Be- 
low the transition layer, z/s ^ 0.2, the incident FUV ra- 
diation is highly attenuated and the chemical timescales 
are controlled by two-body reactions, which are sensitive 
to the temperature and density. The chemical timescale 
in this region is generally longer than the photochemical 
timescale for the interstellar radiation field, 10 3 yr. The 
accretion and mixing timescales become comparable to 
this limit at 0.5 and 10 AU, respectively. The disk wind 
timescales are longer than the chemical timescales for all 
z/s for the entire disk in this model. 

2.7. Infrared molecular line spectra 

Using the results of our chemical evolution calcula- 
tions, we compute emergent mid-infrared spectra from 
the disk. We include the dust continuum, which is the 
dominant source of opacity in protoplanetary disks, along 
with molecular line emission and absorption, assuming 
local thermodynamic equilibrium (LTE). 



Accretion, mixing and wind velocities 



Heinzeller et al. 

Chemical and physical timescales x 
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photochem. 
accretion 
mixing 
wind 



s [AU] 
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Fig. 2. — Velocities (left) and timescales (right) for the transport processes of accretion (radial), mixing (vertical) and wind (vertical). 
Typical timescales of the photochemical reactions are also shown in the right panel (see text for details). Solid lines denote midplane 
values, i.e., z/s = 0, (accretion, mixing, wind), or the timescale of the photochemistry for the interstellar radiation field (photochem.), 
respectively. Dashed lines correspond to z/s = 0.2 and dotted-dashed lines to z/s = 0.5. 



The dust continuum opacities are taken from the 
Jena opacity tables, usin g routines contained within 
the RATRAN package rtOssenkopf fc Hennmgl 119941 : 
IHoge rhciide fc van der Taki l2000fT ~ In accordance with 
the dust properties of the disk model, we use the ta- 
bles for spherical dust grains (parameter "e5" in the 
RATRAN routines) with thin ice coating (parameter 
"thin"). 

The HITRAN'2004 molecular spectroscopic database 
(jRothman et al.ll2005D provides a large set of energy lev- 
els and line transitions in the infrared for 39 molecules, 
focussing on terrestrial atmospheric species. Neverthe- 
less, the overlap of species of interest between the HI- 
TRAN and UMIST databases is sufficient for our pur- 
poses. Since isotopes are not currently included in our 
chemical reaction network, we neglect spectroscopic in- 
formation for minor isotopes in the HITRAN database. 

The major issue concerning the calculation of line in- 
tensities in the infrared is the limited availability of col- 
lisional coefficients, in particular, for the large set of lev- 
els and transitions contained in the HITRAN database 
(e.g., more than 180000 transitions for the main isotope 
of CH4 for wavelengths between 1 and 40 /Win). Since the 
main purpose of this work is the calculation of the disk's 
chemical evolution, a detailed non-LTE radiative transfer 
calculation of the emergent spectra is beyond the scope 
of this paper. We therefore assume local thermal equilib- 
rium to model the line emission and absorption strengths 
for the species of interest. For a species k, the transition 
j —> i (Ei < Ej) is calculated from the standard equa- 
tions 



c (r gas )=Y^9k,i CXP ( ~ |N ) , 



9k i/j 

nfe ^ =nfe Q^)- exp 

jk,ij = hvAkjink,jVk , 



E 



k,i/j 



Kk,ij=hv(B ktij n k ,i - Bk,jiUk,j)Vk 



(22) 
(23) 
(24) 
(25) 



Here, n, j and k denote the level populations, the line 



emissivity and the line opacity, respectively. The Ein- 
stein coefficients A and B arc listed in the HITRAN 
database or determined using the principle of detailed 
balance. The database also provides statistical weights 
g and partition functions Q for a wide range of tempera- 
tures. We calculate Voigt profiles V including natural line 
broadening, pressure broadening and thermal broaden- 
ing. A Doppler profile at a temperature, T gas , is assumed 
for the thermal broadening. Turbulent broadening is ne- 
glected, since it is negligible compared to thermal broad- 
ening for subsonic turbulent velocities. The half-widths 
for natural broadening 7 n and pressure broadening 7 P 
are provided by the HITRAN database at a tempera- 
ture of 296 K. The temperature dependence of pressure 
broadening is approximated by 



7p = 7 P (296K) 



296 K 



gas 



(26) 



where n p is the temperature-dependent exponent in the 
HITRAN database. 

We adopt a simple, one-dimensional linc-of-sight inte- 
gration to calculate the emergent disk spectra for a disk 
viewed face-on, i. e., with an inclination i = 0, relative to 
the observer. The radial resolution for the spectral calcu- 
lation is identical to the chemical model. Along each ray, 
the intensity is integrated following the one-dimensional 
equation of radiative transfer, accounting for the thermal 
dust continuum emission and the molecular line emission, 



3v fcj/Iu 



(27) 



Quantities on the right hand side of (|2"T|) are evaluated 
in the moving frame of a fluid element of the observed 
object, while those on the left hand side are given in the 
rest frame of the observer (denoted with '). The observed 
frequency is given by v' = v (1 + (v/c)cos9), where v 
is the velocity of the fluid element and 9 is the angle 
between the directions of the velocity and the observer. 
The emission and absorption coefficients arc given by 



E E 

k Ei<Ej 



3 k, 



^v.dustB^ (Td us t ) 



(28) 



Chemical evolution of protoplanetary disks 



9 



k Ei<Ej 



(29) 



where B u is the Planck function and K v ^ ust is the specific 
dust absorption. 
The total emission from the disk is calculated as 



(30) 



A/, disk = 87T 2 / sdsl y f(s). 
The full spectrum is obtained from 



v' . disk 



-4:TT 2 R 2 



B v '(T e g) + Ivy', disk • 



(31) 



Here, the stellar emission is modeled as blackbody emis- 
sion with temperature T c g = 4000 K and stellar radius 
R = 2i?Q and is added to the dust continuum and molec- 
ular line emission from the disk. 

For the line-of-sight integration, a volume-weighting 
method is adopted to interpolate the temperature and 
molecular densities before computing the emission and 
absorption coefficients. In doing so, potential effects of 
the relatively coarse grid on the emerging line emission 
are reduced. 

3. RESULTS AND DISCUSSION 

In the following, we discuss the disk chemical evolution 
results and spectral calculations for the different models 
considered in this study. First, we investigate the in- 
fluence of the H2 formation rates on grains on the disk 
chemistry (Section l3Tj) . In Sections l3~2l and f3T3l we study 
the specific effects of viscous accretion, turbulent mixing 
and disk winds. Finally, in Section 13.41 we discuss how 
our results compare with observations of protoplanetary 
disks. 

3.1. H2 grain formation rates 

Recently, iGlassgold et al.l (|2009f ) studied the effect of 
the H2 formation rates on dust grains on the water abun- 
dances at the disk surface of the inner disk. We revisit 
this effect and also study the influence on the infrared 
molecular line spectra. 

In model A, we use th e simp le approximation from 
iGould k Salpeterl (|1963f ) and iHollenbach k Salpeterl 
(|1971l) . In model B, th e adopted H 2 grain form ation 
rate is that according to ICazaux k Tielc ns (2002]) (see 
Section 12.4ft . In both models, the effects of accretion, 
turbulent mixing and disk winds are ignored. The chem- 
ical evolution of the disk is calculated over 10 6 yr for a ra- 
dial range of 0.5-300 AU. For both models, the chemical 
abundances everywhere in the disk have reached steady 
state by approximately 10 5 yr. 

In Figure El we show the resulting gas-phase abun- 
dances for molecular hydrogen, carbon monoxide, wa- 
ter and the hydroxyl radical for the inner 100 AU of the 
disk. In model B, the formation of molecular hydrogen 
is more efficient, particularly in the upper layers of the 
disk. In the warm molecular layer, i.e., the transition 
layer where the temperature rises and the density de- 
creases steeply, H2 is prone to continuous dissociation 
due to thermal reactions and incident radiation, and the 
enhanced grain formation rates in model B are able to 
maintain higher fractional abundances. Similarly, in the 



photodestruction region, i.e., the surface layers of the 
disk, H2 is more abundant by a factor of 10-50 between 
5 and 100 AU. At even smaller radii the strong stellar 
UV radiation destroys H2 on timescales shorter than it is 
formed even in model B. In the lower layers and outer re- 
gions of the disk, the differences between model A and B 
become negligible, since here, H2 accounts for essentially 
all of the hydrogen nuclei. 

The differences in the hydrogen chemistry also affect 
the resulting abundances of predominant molecules in the 
inner disk, particularly CO and H2O. These molecules 
are commonly found in the planet- formin g regions of pro- 
toplanetary disks (c. f. JSalvk et al.|[2008l ). Their modeled 
abundances depend strongly on the amount of molecu- 
lar hydrogen at the disk surface and are therefore very 
different for model A and B. The increase of molecular 
hydrogen by a factor of 10-50 between 5 and 100 AU 
in the upper layers translates to an increase in CO and 
H 2 by a factor of ps 1000. Such a strong co rrelation is 
expected theoretically (jGlassgold et al.ll2009() : high frac- 
tional abundances of H2 of ~ 10~ 3 in the surface layers 
produce significant amounts of CO, H2O and OH through 
a sequence of neutral radical reactions: 



+ H 2 
OH + H 2 
C + OH 



H + OH 
H 2 + H 
CO + H. 



IGlassgold et al.1 (|2009f ) show that the ratio z(H 2 0)/x(0) 
depends quadratically on the ratio x(H 2 )/a;(H + ). Our re- 
sults confirm this correlation, since the O and H + abun- 
dances are practically identical between model A and B. 

The molecular column densities are only slightly 
changed compared with the fractional abundances, since 
the the latter are altered only at the disk surface where 
the density is low (Figure [4|). Independent of the H 2 for- 
mation rate, gaseous CO traces the total disk mass, with 
a constant ratio A r co/A r H 2 = 10~ 4 at radii < 100 AU 
where the dust temperature is high enough so that CO 
does not freeze out. Gas-phase water is fairly abundant 
at radii < 3AU (i.e., within the snow line of the disk) 
with Nn 2 o/Nco = 0.25. The column density of water 
is increased by about one order of magnitude between 5 
and 50 AU, which is due to the large overabundance in 
the upper layers of model B. 

The two H 2 formation rate models affect the abun- 
dances of the observable molecules in the warm molec- 
ular layer and the photodestruction layer, from where 
infrared molecular line emission originates. In Figure [5J 
we display the emergent spectra for a disk seen face on 
(inclination i = 0). The spectral resolution is set to 600 
to match the Spitzcr/IRS specifications. The molecular 
line emission is calculated for the following molecules be- 
tween 1 and 40/xm: CO, H 2 0, OH, C0 2 , CH 4 , NO, S0 2 , 
NH 3 , C 2 H 2 , N 2 , H 2 S, and NO+. The spectra differ sub- 
stantially between the two models. Model A produces 
strong and narrow emission lines of OH between 10 and 
40/xm, due to high abundances in the surface layers of the 
disk. On the contrary, H2O shows broad emission bands 
between 5 and 10/xm, but only weak and narrow emission 
and absorption features at larger wavel engths, in con- 
trast to the observations (c. f., Figure 1 in lCarr k Naiital 
12001 . The CO emission band at 4-5 /tm is present, as 
well as an absorption dip of CO2 at 15/xm. The absorp- 
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Fig. 3. — Gas-phase fractional abundances of H2, CO, H2O and OH for models A and B 
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tion signatures are generated at radii ;$ 2 AU, where the 
temperature increases near the midplane due to viscous 
heating (c. f., the inlay in Figure [1]). Since the abun- 
dances of CO2 and H2O in model A arc high near the 
disk midplane, but low in the disk surface where the gas 
temperature increases towards an observer, absorption 
lines rather than emission lines are generated. Emission 
or absorption from other molecules included in the spec- 
tral calculation (e.g., NH3, CH4) is very weak in the 
spectrum. 

Model B, on the other hand, produces a spectrum with 
enormous emission features from OH, H2O and CO. The 
CO2 absorption feature is not affected by the H2 forma- 
tion model, since it is destroyed very efficiently in the 
upper layers of the disk, where the main differences in 
the molecular abundances arise between model A and 
B. Comparing the two-dimensional abundances of OH, 
H 2 and CO in model A and B (Figure [3]) allows us 
to constrain the origin of the emission lines to roughly 
3-20 AU, where the gas temperature is high enough to 
contribute to the total flux (c. f., Figure [J). Such strong 
emission features are certainly in disagreement with ob- 
serv ations. The enhanced form ation rate of H2 on grains 
from lCazaux fc Tielensl ((2002f ) might simply produce too 
much H2 due to the choice of the variety of parameters 
used for this model, however, we should note here the 
limitations of our simple dust model. Firstly, we assume 
that gas and dust are well mixed in the disk, with a 
constant mass density ratio of p ga s/ Pdust ~ 100. Numer- 
ous studies show that dust grains will settle towards the 
mid plane instead of residi ng in the upper disk layers (see, 
e. g.. lDominik et al.ll2007l for a review). Less dust grains 
in these layers lead to lower H 2 abun dances and thus to 
lower abundances of water and OH. lAikawa &: Nomural 
(|2006| ) demonstrated that less dust grains will also lead 
to a lower gas temperature at the disk surface, which 
reduces the line intensities. Since the same line of ar- 
gument applies to model A, it implies even lower abun- 
dances and thus weaker emission lines for H2O in this 
model. Secondly, the spectral calculations in our study 
assume LTE level populations, w hich may overestimate 
the strength of the emission lines. iMeijerink et al.l ()2009l ) 
studied the effects of non-LTE radiative transfer on the 
molecular emission lines of water: they assume fractional 
H2O abundances between 10 -10 and 10~ 4 throughout 
the disk. Interestingly, they conclude that a strong de- 
pletion at s < 1 AU is required to account for observed 
emission line spectra, and such a depletion is found in 
our model A. There, the H2O abundance is low in the 
surface layer, because of the efficient destruction of wa- 
ter due to the stellar irradiation and th e high tempera- 
ture. Most important, however, is that iMciicrin k et al.1 
(|2009T l find drastic differences between LTE and non-LTE 
spectral calculations: the LTE calculations overestimate 
the strength of the emission lines by a factor of 10 for 
the strongest transitions (i.e., the largest Einstein coef- 
ficients). In model B, the line emission is generated in 
the upper layers of the disk, where the gas number den- 
sity is of the order of 10 8 cm~ 3 . The water line transi- 
tions betwee n 10 and 40/xm report ed from Spitzer obser- 
vations (e. g-. ICarr k, Na iita 200§) are mainly rotational 
transitions in the ground vibrational state and in the first 
vibrati onal state with critic al densities of 10 12 cm~ 3 or 
larger ([Meiierink et al.l l2009) . Vibrational transitions be- 



tween the two states that fall into this wavelength range 
have even higher critical densities. Hence, the conditions 
for LTE are partly violated and a proper non-LTE cal- 
culation of the disk spectra should be considered. In the 
following sections, we ignore the enhanced H2 formation 
rates of model B and study the effects of accretion, tur- 
bulent mixing and wind on model A only. 

3.2. Accretion, mixing, disk winds - single column 
solutions 

To investigate the effects of including physicochemical 
interaction, we first focus on one vertical column of the 
disk at radius s = 1.3 AU. Here, the temperatures of gas 
and dust are high enough, in particular in the upper lay- 
ers of the disk, to allow for an active chemistry. Further, 
near- and mid-infrared observations derive characteris- 
tic radii for molecular line emission around 1 AU (c. f., 
Section [374]). 

In the following, we discuss in detail the modeled 
molecular abundances for a number of species of interest. 
For a better understanding, in particular for interpreting 
the resulting fractional abundances xu = nk/n to t, we 
display the physical properties of the disk column at this 
radius in Figure O The gas and dust temperatures are 
tightly coupled in the dense midplane with values around 
200 K. Hence, practically all species have desorbed from 
the dust grains and the ice abundances are negligible. In 
the warm molecular layer {z w 0.2 AU), the gas and dust 
temperatures begin to decouple due to the lower den- 
sity. In the photodestruction layer, the gas is heated up 
to 4000 K because of the strong stellar irradiation, while 
the dust temperature is around a factor of 10 lower. 

In Figure we plot the fractional abundances as a 
function of vertical distance from the midplane for the 
model without physicochemical interaction (model A) 
and for models with either accretion, turbulent mixing 
or disk winds (models ACR, MIX, WND). The gas-phase 
abundances are shown together with the - mostly neg- 
ligible - ice abundances. All four models make use of 
the simple H2 formation rates on grains (Eq. ([2])), while 
model B (not displayed) includes the enhanced H2 for- 
mation rates on grains (Eq. ((4|)) without physicochemical 
interaction. Table [3] at the end of this section compares 
the column densities of the most important species at 
radius 1.3 AU for all models A, B, ACR, MIX, WND. In 
the following, we discuss in detail the effects of turbulent 
mixing on the vertical chemical stratification at a radius 
,s = 1.3 AU, since this model has the largest effects on 
the chemical structure of the disk. 

3.2.1. Turbulent mixing 

Turbulent mixing has important effects on the disk's 
chemical composition, in particular in the upper layers. 
Generally, diffusive mixing tries to homogenize the frac- 
tional abundances of molecules over the depth of the disk. 
Figure [7| illustrates this effect, for example in the case of 
H2, the fractional abundances are slightly increased in 
the transition region, with only negligible changes in the 
column density (c. f., Table [3] at the end of this section). 

This effect is more pronounced for CO2: in model A, 
CO2 is abundant close to the midplane only, since it is 
removed from the upper layers by reactions with H and 
H + and by photo-dissociation due to the incident UV 
and X-ray radiation: 



12 



Model A: column number densities 



Heinzeller et al. 

Model B: column number densities 
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Fig. 4. — Column densities for the species of interest in models A and B. Solid lines correspond to gas-phase species, dashed lines to ices. 



Model A: disk spectrum for inclination ;'=0 



Model B: disk spectrum for inclination !=0 
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Fig. 5. — Emergent LTE disk spectra for models A and B for dust continuum and molecular line emission. 
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Fig. 6. — Gas and dust temperature T gas and T,j us t, and total 
number density ntot as function of vertical distance from the disk 
midplane at a radius s = 1.3 AU. 



MIX, CO2 diffuses upwards and its abundance is reduced 
near the midplane. 

The H2O, OH, and CO abundances are correlated with 
that of molecular hydrogen and are enhanced in model 
MIX in the region where H2 is more abundant. The ef- 
fect of diffusion of H 2 leads to reduced fractional abun- 
dances for z < 0.15 AU and increased abundances be- 
tween 0.15 AU and 0.2 AU. Water is destroyed easily 
at high temperatures and through photodissociation by 
UV and X-ray photons. The H2O abundance has a dip 
at z w 0.22 AU, where both the gas temperature and 
density are relatively low (c. f.. lWoitke et alJr2009bl Fig- 
ure 1). Here, UV and X-ray photodestruction become 
strong, while at the same time the temperature is still too 
low to activate the main formation reactions for H2O, 



H 2 + 
OH + H2 



OH 4 
H 2 



H 

-H. 



H- 
H+ 
C0 2 



C0 2 
C0 2 

- 7uv 



C0 2 + 7X-ra 



CO + OH 
HCO+ + O 
CO + O 
CO + O. 



At z w 0.2 AU, the reactions with H + and the X-ray 
photochemistry are dominant, while in the upper layers, 
the UV photochemistry is most important. In model 



In the case of methane, diffusion also increases the frac- 
tional abundances in the surface layers. At the same 
time, the formation of CH4 close to the midplane is effi- 
cient enough to maintain its initial fractional abundance, 
leading to an overall increase in the column density of 
methane. For z > 0.3 AU, methane is destroyed quickly 
by UV and X-ray photons. 

A similar mechanism is responsible for the considerable 
increase in ammonia for z < 0.3 AU. The formation of 
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z [AU] z [AU] 

Fig. 7. — Fractional abundance profiles of gas-phase species (solid) and ices (dashed) as a function of vertical distance from the midplanc 
at s = 1.3 AU in model A (black), and in models with accretion (ACR, blue), mixing (MIX, red) or disk winds (WND, green). 
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Fig. 7. — (cont.) Gas and ice fract. abundances at s = 1.3 AU. 



ammonia near the disk midplanc is driven by a change 
in the oxygen chemistry, with O2 being overabundant in 
model MIX (see also Section [3~3)h The differences in the 
oxygen abundance between model A and model MIX is 
caused by the continuous diffusion of O2 from the warm 
molecular layer towards the midplanc. This implies a de- 
crease in the column density of HCN, whose formation 
is competing with NH3 for the available nitrogen. Since 
the number of nitrogen atoms over the full column is con- 
served and since the formation of NH3 is more efficient, 
HCN is depleted by one order of magnitude close to the 
midplane, compared with that in model A. Its column 
density decreases by approximately the amount required 
for the formation of the additional ammonia molecules. 

The most pronounced effect arises in the case of 
methanol as a consequence of the increase in ammonia. 
Its fractional abundance increases for z < 0.3 AU by one 
to six orders of magnitude, resulting in a column density 
of 1.2 • 10 19 cm -2 compared with 4.7 • 10 17 cm~ 2 in model 
A. This is du e to an enhanced productio n of CH3OH by 
the reaction (jRodgers fc Charnlevll200ll ) 



CH 3 OH++NH 3 



NH4 + CH3OH . 



In Figure [3 we display the ammonia and methanol abun- 
dances for an additional model calculation with diffusive 
mixing, but where this reaction is switched off. The am- 
monia abundances remain unchanged, but the methanol 
abundances drop by one order of magnitude, which is 
comparable to the midplanc abundances in model A. 

Significant changes are also found for acetylene, C2H2. 
At z = 0.2 AU, it is mainly produced from NH3 by the 
reactions 



NTU 



NH+, 



NH3 + C2H+ — > C 2 H 2 

NH 3 + C 2 H+ — > C 2 H 2 

as well as by a reaction of C and CH3 at z < 0.2 AU, 

C + CH3 — > C 2 H 2 + H . 

Switching off these reactions leads to a significant de- 
crease in C2H 2 (c. f., Figure [7|). Hence, acetylene profits 
from the increase in ammonia abundances and is over- 
abundant by two orders of magnitude at z < 0.2 AU. 
Higher up in the disk, it is destroyed efficiently by photo- 
chemical reactions and reactions with atomic hydrogen. 

Finally, we inspect the sulphur-bearing species SO2 
and H 2 S. Initially, sulphur is locked in H 2 S only, with 
a fractional abundance of 10~ 6 (c. f., Table [2|). H2S is 
easily destroyed by atomic hydrogen due to the small ac- 
tivation barrier of this reaction. Provided the O2 and 
OH abundan ces are sufficient l y high, this drives a chain 
of reactions (jCharnlevi Il997t (Millar et all 11997( 1 which 
forces large amounts of H2S into SO2 in both models 
at z < 0.25 AU. The now highly abundant SO2 near the 
midplane is destroyed on timescales longer than its for- 
mation timescales and therefore provides a steady reser- 
voir to maintain the enhanced H2S abundances in com- 
parison to model A. At the midplane, the destruction of 
H 2 S is less effici ent, but still sufficien t to produce OCS 
from CO and S (jHatchell et al.lll99^ 1. After 10 6 yr, sul- 
phur is locked in SO2 and OCS in equal amounts. 

Contrarily, in model A, the formation of SO2 is ineffi- 
cient in the lower layers due to the missing OH and the 
missing O2 at z 0. Thus, the slower reaction path 
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of S with CO is preferred and OCS is built up over the 
lifetime of the disk, eventually locking all sulphur. 

3.2.2. Disk winds 

The resulting chemical abundances for model WND 
are displayed in Figure [7] For all species, the effects of 
the disk wind on the vertical stratification of the chem- 
ical species are much less pronounced than in the case 
of turbulent mixing. In the cases of CO2, NH 3 , CH 3 OH 
and C2H2, slight increases in the fractional abundances 
are found around z = 0.2 AU, in particular where narrow 
dips occur in the abundances from model A. Protected 
by the upper layers against the UV and X-ray irradia- 
tion, molecules transported upwards from the midplane 
fill these sinks on timescales comparable to their removal 
by chemical reactions. 

The reason for the inefficiency of the disk wind lies in 
the steady disk model assumption considered here: to 
avoid disk dispersal over the lifetime of the disk, which 
would contradict the steady state assumption, we in- 
troduced a sc aling factor n = 1 0~ 3 f or the resulting 
wind speeds of lSuzuki fc Inut suka (20091) at the midplane 
layer. Hence, the upward velocities near the midplane are 
low in model WND, compared to the turbulent velocities 
in model MIX (c.f., Figure 0). At radii s < 10 AU, the 
wind velocities become very large in the photodestruc- 
tion region of the disk (z/s > 0.3), due to the continuity 
equation in vertical direction in steady state. But there, 
the timescales of the photochemistry are extremely short 
so that the effect of the disk wind is small (c. f., Figure[5]). 

For a scaling factor h = 10~ 3 , the midplane grid cell 
at radius s = 1.3 AU evaporates in about 2 • 10 5 yr, cor- 
responding to a mass loss of 20% of the total mass con- 
tained in this disk ring. Hence, in a stationary model, 
x should be lower than 10 -3 , which would lead to even 
less efficient disk winds. 

We thus conclude that in a steady disk model, the disk 
winds have only negligible effects. Due to the neces- 
sary downscaling of the wind speeds by a factor of 1000 
[k = 10 -3 ), the vertical mass transport becomes small 
and causes only small changes in the disk's chemical com- 
position. Nonetheless, it is worth noting that in a time- 
dependent physical disk model, where evaporation and 
dynamical refueling by the accretion flow are allowed, 
the disk wind may become significant. The simple model 
considered here suggests that it will affect molecules that 
also show an enhancement in the upper layers in the dif- 
fusive mixing model (e.g., NH3, CH3OH). 

3.2.3. Viscous accretion 

The chemical abundances from model ACR at a radius 
of 1.3 AU are displayed in Figure [7] As in the previous 
cases, at the disk surface photochemistry dominates and 
the accretion flow does not affect the chemical composi- 
tion. Close to the disk midplane, the accretion time scale 
is comparable to the chemical timescale and the effects 
of the radial inflow of material become visible. Whether 
the accretion flow leads to an increase or a decrease in 
the abundance of a particular species depends on the 
formation and destruction timescales of that species. In 
the case of H2O, for example, the chemical reactions are 
fast enough to compensate for the effects of the accretion 
flow and the abundances and column density are similar 
to those from model A. 



TABLE 3 

Calculated column densities JV[cm -2 ] at radius 

! = 1.3 AU OF SELECTED SPECIES FOR THE DISK MODELS 
CONSIDERED IN THIS STUDY. 



species 


A 


B 


ACR 


MIX 


WND 


H 2 


1.7e+25 


1.7e+25 


1.7e+25 


1.7e+25 


1.7e+25 


CO2 


2.5e+19 


2.6e+19 


2.5e+19 


8.3e+18 


2.4e+19 


H2O 


4.2e+20 


4.2e+20 


4.0e+20 


3.5e+20 


4.2e+20 


OH 


1.5e+15 


9.3e+15 


1.5e+15 


5.3e+15 


1.8e+15 


CO 


1.6e+21 


1.6e+21 


1.6e+21 


1.7e+21 


1.6e+21 


CH 4 


5.4e+18 


5.4e+18 


6.5e+17 


r.ic+18 


5.7c+18 


NH 3 


2.5e+20 


2.5e+20 


3.4e+20 


3.0e+20 


2.6e+20 


HCN 


4.3e+19 


3.8e+19 


l.lc+19 


4.7e+18 


4.0e+19 


CH3OH 


4.7e+17 


4.8e+17 


1.4e+16 


1.2e+19 


5.6e+17 


C2H2 


1.6e+16 


1.6e+16 


1.3e+14 


2.2e+18 


1.6e+16 


SO2 


2.7e+17 


2.6e+17 


1.2e+19 


1.4e+19 


2.9e+17 


H 2 S 


9.6e+14 


8.3e+17 


6.2e+15 


2.8e+18 


9.4e+14 


OCS 


3.3c+19 


3.2e+19 


2.0e+19 


1.6c+19 


3.3c+19 



On the contrary, NH3 shows higher abundances at z < 
0.15AU. As discussed for the case of diffusive mixing 
above, the total amount of nitrogen atoms is conserved 
and determined by the initial abundances of N2 and NH3 
(c. f., Table[2]). An increase in the column density of NH3 
from 2.5 • 10 20 cm -2 in model A to 3.4 • 10 20 cm~ 2 in model 
ACR leads to a decrease in HCN from 4.3 • 10 19 cm~ 2 to 
1.1 • 10 19 cm~ 2 at this radius. 

In contrast to the results from model MIX, methanol 
shows lower abundances close to the disk midplane, de- 
spite the increase in ammonia. Since its fractional abun- 
dance is even lower than those from model A, the reason, 
therefore, is not only the lack of vertical diffusion, but the 
fact that the formation and evaporation timescales (from 
dust grains) of methanol at around 10 AU are longer than 
its destruction timescale (c.f., Section l3~3|) . The accre- 
tion flow propagates this depletion inwards over time. 
Compared to model A, this reduces the CH3OH column 
density from 4.7 • 10 17 cnT 2 to 1.4 • 10 16 cm" 2 . 

One should bear in mind that the location of the 
snow line for methanol and hence its gas-phase abun- 
dance is sensitive to its binding energy on dust grains, 
which varies across the li terature. Here, we adop t 
the theoretical va lue from Hascgawa & Hcr bst] (| 1993T) . 
£ bind = 2140 K. iSandford fc Allamondolal (119931) and 
more recently iBrown fc BolinaT i 20071) find much higher 
values between 4000 and 5000 K based on the results 
of temperature-programmed desorption (TPD) experi- 
ments on pure methanol ice. Higher binding energies 
would imply lower gas-phase abundances and a shift of 
the snow line to smaller radii. Here, the NH3 abundances 
are higher and the formation of methanol more efficient, 
which would imply higher gas-phase abundances within 
the snow line. 

A similar decrease in the abundances in the innermost 
disk can be found particularly for C2H2, where the ac- 
cretion flow removes a considerable fraction of acetylene 
close to the midplane (c. f., Table [3] for the resulting col- 
umn densities of C2H2 and other species). 
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3.3. Accretion, mixing - full disk solutions 

In the previous section, we inspected in detail the 
chemical abundances for the different disk models as a 
function of height at a single radial position. The chemi- 
cal composition varies between the individual models due 
to mass transport in the radial and vertical directions. 
The results thereby depend strongly on the species con- 
sidered, in particular, on its formation and destruction 
timcscalcs relative to the dynamical timcscalcs. Since 
these timescales vary with radial distance from the star 
for accretion, turbulent mixing and disk winds, we expect 
significant changes in the effects of these mass transport 
phenomena at different radial positions in the disk. 

Here, we focus on the radial accretion and vertical mix- 
ing and neglect the disk wind. At a radius, s = 1.3 AU, 
the influence of the disk wind on the chemistry is small, 
compared with diffusive mixing, due to the low wind 
speed required to impede the evaporation of the disk 
in our model. At larger radii, the wind velocities are 
even smaller, leading to longer timescales of the upward 
wind motion. At the same time, however, the chemical 
timescales increase (c.f., Section [2.6. 4|) . so that we do 
not find significant changes to the disk chemistry due to 
the inclusion of a disk wind. 

We compare the resulting fractional gas-phase abun- 
dances and the gas and ice column densities for model 
ACR and model MIX with model A in Figure U The 
results for water differ only slightly between the mod- 
els, and the only difference arises in a narrow stripe 
z/s S3 0.2-0.4, where the gas-phase abundances of wa- 
ter are increased in model MIX due to the enhanced H2 
abundances (see also the discussion in Section l3~Tj) . This 
leads to an increase in the gas-phase column density be- 
tween 2 and 10 AU. Since the H2O abundances are in- 
creased in the transition layer, turbulent mixing affects 
the strength of the H2O line emission from the disk. This 
can be seen in Figure [9l where we compare the emerging 
disk spectra for model A, model ACR and model MIX. 
Apart from H2O, the species dominating the molecular 
line emission are CO and OH. Both show slightly stronger 
emission lines in model MIX and in model ACR, although 
the effect is more pronounced in the case of diffusive mix- 
ing. 

In the case of CO2, turbulent mixing reduces the gas- 
phase column density inside 10 AU: turbulent mixing 
leads to the diffusion of CO2 from the lower, protected 
layers into the upper disk, where it is efficiently destroyed 
by incident radiation, compared with its production. For 
radii < 1 AU, CO2 is strongly depleted. Viscous accre- 
tion, on the other hand, enhances the CO2 gas-phase 
abundances at radii < 1 AU due to a net inward transport 
from radii < 10 AU. This leads to changes in the disk 
spectra for both models ACR and MIX: inside 2 AU, the 
gas and dust temperature first decrease with increasing 
distance from the midplane (see Figure [T|). This is due 
to the viscous heating near the midplane. Since CO2 is 
concentrated near the midplane and since it is abundant 
in model ACR (depleted in model MIX), the absorption 
feature at 15/mi is enhanced (reduced). 

Ammonia shows the largest variations between the dif- 
ferent models. Compared to model A, model ACR and 
model MIX lead to a similar increase in the gas and ice 
abundances between 1 and 10 AU. The increase in am- 



monia is a result of a fundamental change in the oxy- 
gen chemistry. In model ACR and model MIX, atomic 
oxygen is reduced by orders of magnitude between 2 and 
10 AU, with an increase in O2 and a decrease in OH. The 
drastic decrease in atomic oxygen practically switches off 
the the main destruction reaction of NH3 , 

O + NH3 — > OH + NH 2 . 

The overabundant regions are located below the warm 
molecular layer in model ACR and reach slightly into it 
in model MIX (see also Figure [7J . Hence, ammonia is 
found in absorption in the spectrum of model MIX due 
to the same reasons as CO2 in model ACR. 

Another species of interest is methanol, whose gas- 
phase abundance is high in regions where the ammonia 
abundances are high in all the models. In model MIX, 
this implies an increase in methanol within 5 AU com- 
pared to model A. Viscous accretion increases its abun- 
dance between 2 and 5 AU and decreases it further inside. 
As in model MIX, CH3OH is produced efficiently from 
NH3, however, its gas-phase formation and evaporation 
from dust grains at 5-10 AU is slower than its destruc- 
tion and the accretion flow transports less material into 
the inner disk than is removed from it. Methanol is lo- 
cated below the transition region of the disk, even in 
model MIX, hence, the line opacities are so small that 
the absorption features are negligible in the spectra. 

The abundance profile of acetylene can be explained 
as follows: in model MIX, its abundance is increased be- 
tween 1 and 10 AU. Further inside, the diffusion leads to 
a depletion of C2H2. In model ACR, C2H2 is increased 
between 2 and 10 AU and decreased further inside, simi- 
lar to CO2. Since the overall abundance of C2H2 is low 
and since it is located below the transition region, the 
absorption lines are negligibly small in the spectra of the 
disks in model A and model ACR. The stark increase in 
C2H2 in model MIX leads to a weak absorption feature 
seen at 13-13. 5/xm (Figure [9}. 

3.4. Theory versus observations 

We refer to the observations of t he protoplanetary disks 
of AATauri (iCarr fc Naiital l2008l). DRTauri and AS 205 
(jSalvk et all 12008ft . AATauri is a classical TTauri star 
with mass .8M(T), accretion rate 1 0~ 8 M ( ^yr~ 1 and incli - 
nation 75° (jBouvier et al.lll~999l ). ICarr fc Naiital pOOl ) 
observed the disk around AA Tauri with the Infrared 
Spectrograph (IRS) on the Spitzcr Space Telescope over 
the wavelength range of 9.9-37.2/im. They found a rich 
spectrum of molecular emission lines dominated by ro- 
tational transitions of H2O and OH. Between 10 and 
15/iin, they detected ro- vibrational emission bands of 
C2H2, HCN and CO2, plus the atomic [Nell] transition 
at 12.8/im. Fundamental ro-vibrational emission bands 
of CO have been detected in Keck/NIRSPEC observa- 
tions betwee n 3 and 5^m. Fitting model spectra to the 
observations. ICarr fc Naiital (|2008| ) derive a high column 
density for these species and characteristic emission radii 
in the inner disk (s < 2.5 AU). 

DRTauri and AS 205, both classical TTauri stars, 
have been observe d with IRS on S pitzer and with 
Keck/NIRSPEC by [SalyFeFin (llOOl). DRTauri has 
a mass of 0.76M Q and a variable accretion rate of [0.3- 
79] • lO -7 M0yr -1 , much higher than in our model or that 
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Fig. 8. — Two-dimensional fractional gas-phase abundances and column densities (solid: gas; dashed: ice) of H2O and CO2 for model A 
(black), model ACR (blue) and model MIX (red). 
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Fig. 8. — (cont.) Two-dimensional fractional gas-phase abundances and column densities (solid: gas; dashed: ice) of NH3 and CH3OH 
for model A (black), model ACR (blue) and model MIX (red). 
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Fig. 8. — (cont.) Two-dimensional fractional gas-phase abundances and column densities (solid: gas; dashed: ice) of C2H2 for model A 
(black), model ACR (blue) and model MIX (red). 



observed in AA Tauri. AS 205 is the primary component 
of a triple system with mass 1.2M© and accretion rate 
7.2 ■ lO" 7 M yr- 1 . The inclinations are 67° for DR Tauri 
and 47° for AS 205A, respectively. Both systems show 
a large number of transitional emission lines from water 
and require high temperatures and column densities to 
explain this emission, similar to AA Tauri. Further, CO, 
OH and CO2 are required to fit the observed spectra in 
the near- to mid-infrared. In the following, we discuss 
how our models relate to, and compare with, these ob- 
servational results. 

The calculated spectra for our models show emission 
and absorption features different from the observed spec- 
tra (Figure [9]). Firstly, the molecular line emission in 
our models is dominated by rotational transitions of OH. 
Weak emission from H2O is detected in model MIX. In 
model ACR and in model A, water is found mainly in ab- 
sorption. Model WND leads to results similar to model 
A due to the limitations of our disk model. These mod- 
els make use of the simple H2 formation rates on grains, 
Eq. ([2]). Model B, which accounts for the enhanced H2 
formation rates on grains, Eq. (|4|), leads to huge emis- 
sion lines from OH and H 2 0. The near- infrared emission 
bands of CO at 5/xm are reproduced in model A, model 
MIX and also in model ACR. 

The modeled line intensities depend strongly on the gas 
temperature and molecular abundances profiles, i. e., the 
detailed thermal balance and the chemistry. It is there- 
fore more favorable to compare the molecular column 
densities from our model calculations to the observed 
ones. In the inner disk, the observed column densities 



and abundance ratios are derived from the photodestruc- 
tion layer and the warm molecular layer, i.e., from the 
disk surface down to an optical depth of about unity. 
The total molecular column densities in our model are 
mainly determined by the dense midplane layers and are 
therefore larger by orders of magnitude and not compa- 
rable to the observed column densities. Thus, we cal- 
culate the column densities of the observed species from 
the disk surface down to a continuum optical depth of 
unity at each wavelength and compare them to the ob- 
servations (Figure [T0|) . At s = 1.3 AU and wavelengths 
A = {5, 10, 20, 40} /im, the disk becomes optically thick 
at z/s = {0.14,0.15,0.14,0.11}, respectively. For the 
unusually strong OH lines, a continuum optical depth of 
unity corresponds to a maximum optical depth of r ~ 20 
(see discussion below for large discrepancy between the 
models and the observations of OH lines). For all other 
molecules, the maximum optical depth of the lines is 
t - 2. 

If we compare the H2O and CO column densities of the 
upper layers at disk radius 1.3 AU, we find an abundance 
ratio of 0.02 for model A and model AC R, while model 
MIX shows an abundance ratio of 0.18. iCarr fc Naiital 
(12001 der ive an a bunda nce ratio of H2O to CO of 1.3, 
while ISalvk et all ([2008) find water to be less abundant 
than carbon monoxide with an abundanc e ratio for H2O 
to CO of 0.1 for both systems. Note that ICarr &: Naiital 
(|2008| ) calculate the abundance ratio for different radii of 
the H 2 Q and CO colum n densities (2.1 AU and 0.7 AU), 
while ISalvk et al.1 (|2008l ) obtain their results for the same 
radius (3AU). It is important to note that the ob- 
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Fig. 9. — LTE line spectra for model A, model ACR and model MIX for an observer at inclination i = 0. For model MIX, we also 
show the spectrum normalized by the continuum flux in the lower right panel (for a better illustration, the color coding in the normalized 
spectrum is different from the other spectra and listed in lower right legend). 



served systems differ from our model in their central ob- 
jects, disk masses, inclination angles and accretion rates, 
with AATauri being the cl osest match. Further , the 
derived column d ensities in iCarr fc Naiital (|2008| ) and 
iSalvk et al.l (j2008| ) contain considerable uncertainties in 
their values and their characteristic radii. 

The resulting column densities of the upper disk lay- 
ers are displayed in Figure [TO] for model A and for model 
MIX only, since it leads to the largest changes in the ver- 
tical distribution of the chemical species. In model MIX, 
the H 2 to CO abundance ratio is of the order of 0.2 
inside the snow line, while in model A, it is reduced by 
at least a factor of 10. Model MIX reproduces the H 2 
abundance in AATauri at a radius of 2.1 AU. In the 
case of DRTauri and AS 205, the central star is more lu- 
minous than in AA Tauri and in our model, hence, the 
snow line is located further outwards and gas-phase H 2 
remains abundant with derived column densities of the 
same order as in our model or in AA Tauri at radii of 
2 AU. The modeled CO abundances agree with the ob- 
servations for both models. Both models fail to repro- 
duce the high OH abundances, although the calculated 
line emission is much stronger than that of the observa- 
tions. This adds further support to the call for non-LTE 
spectral calculations and gas temperature profile calcula- 
tions accounting for dust-grain settling and grain growth 
(c. f., Section HO]) . Since the binding energy of OH on 
dust grains is low and considering it is robust against 



UV radiation, efficient turbulent mixing can increase the 
abundances in the upper layers of the disk (c. f., Fig- 
ure [TO]). 

C0 2 is found in absorption instead of in emission in 
model A, model ACR and model MIX. The absorption 
is strongest in model ACR and weakest in model MIX. 
Contrary to this, the observations detect C0 2 in emission 
at a characteristic radius of 1.2 AU in AATauri. It is 
also detected in emission towards DR Tauri and AS 205 
(with large uncertainties on the column densities and the 
emission radii). The column densities derived from the 
observations of AA Tauri are lower than in model A and 
model MIX, although the latter one gives a better fit. 
Since C0 2 is destroyed quickly in the upper layers of the 
disk, stronger mixing will lead to lower C0 2 abundances. 

HCN is detected in weak emission in AATauri 
and a number of oth er sources (see, for example, 
IPontoppidan et all 120101 ) . but not detected in our mod- 
eled spectra. The column density in AA Tauri is slightly 
higher than in our models, and model MIX again gives a 
better fit. Turbulent mixing increases the column density 
in the upper layers, although the total column density is 
reduced. Stronger mixing will help to improve the fit of 
the model calculations to the observations. 

In the case of C 2 H 2 , the characteristic radius in 
AA Tauri could not be determined and was set to that 
derived for HCN. Assuming that the C 2 H 2 emission is lo- 
cated at larger radii (e. g., 1 AU instead of 0.6 AU), model 
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Fig. 10. — Column densities of the gas-phase species in the upper 
layers in model A (dashed) and model MIX (dotted). The derived 
column densities for AA Ta uri, DR Tauri and AS 205 are indicate d 
at their characteristic radii (Carr & Naiita 2008; Salvk et al. 2008). 

MIX fits the observed column density, while model A 
predicts much lower values. Adversely, the spectrum of 
model MIX shows C2H2 in weak absorption instead of 
emission. From the previous discussion in Sections 13.21 
and 13.31 we can assume that stronger turbulence will 
transport C2H2 higher up in the disk, leading to emis- 
sion instead of absorption signatures. 

The abundances of NH 3 arc enhanced in model MIX, 
leading to absorption signatures in the spectrum between 
8 and 13/im. To date, there has been no published detec- 
tion of ammonia from protoplanetary disks, but NH3 has 
been found in absorption towards the hot cor e associated 
with the massive protostar NGC 7538 IRS 1 (jKnez et al.1 
[2009I ). using the Texas Echelon Cross Echelle Spectro- 
graph (TEXES). Observing NH3 emission from disks 
with high resolution and high sensitivity will allow us 
to draw important conclusions on the efficiency of tur- 
bulent mixing in protoplanetary disks. 

Taking into account the model uncertainties, e. g. the 
H2 formation rates on grains, the gas temperature pro- 
files (which are related to the dust model adopted, e.g., 
to the dust growth and settling), as well as the uncertain- 
ties in deriving molecular column densities and charac- 
teristic radii from the observations, further investigations 
are needed to strengthen the significance of our results. 



In addition, the observed objects have inclination an- 
gles of 67° for DR Tauri and 47° for AS 205 A, while our 
theoretical spectra are calculated for a disk seen face-on. 
This affects the observed line emission and the derived 
column densities, since the chemical composition is differ- 
ent along the line of sight. Robust conclusions therefore 
require spectral calculations accounting for the different 
inclination angles. 

4. CONCLUSION 

We studied the chemical evolution of a steady proto- 
planetary disk, investigating different models for the for- 
mation rate of H2 on dust grains and the effects of includ- 
ing viscous accretion, turbulent mixing and disk winds. 
Using our chemical results, we calculated the emergent 
near- and mid-infrared spectrum from the disk, includ- 
ing contributions from the dust continuum and molecular 
line emission. 

Our results show that the abundances of water and 
the hydroxyl radical are very sensitive to the H2 for- 
mation rates on grains, with the resulting spectra dif- 
fering considerably. The enhanced H2 formation rates of 
ICazaux fc Tielc ns ( 2002]) lead to a drastic increase in the 
gas-phase abundances of water and OH in the upper disk 
layers and to strong emission lines. We concluded that 
dust grain settlement towards the midplanc and non-LTE 
calculations of the molecular line emission are required 
to obtain meaningful results for models using these en- 
hanced H2 formation rates. 

We investigated the effects of physical mass transport 
phenomena in the radial direction by viscous accretion 
and in the vertical direction by diffusive turbulent mix- 
ing and disk winds. Due to the steady state assump- 
tion of the underlying protoplanetary disk model, the 
wind velocities are limited to impede the evaporation of 
the disk over its lifetime and consequently, the effects 
of the disk winds on the chemistry are small. On the 
contrary, vertical diffusive mixing has a considerable im- 
pact on the disk chemistry and increases the abundance 
of many observed species in the warm molecular layer 
of the disk. Species which are produced efficiently in a 
particular layer and which are stable in lower (or up- 
per) layers can experience a significant enhancement due 
to the diffusive motion. Including the radial accretion 
flow in the chemical network calculation leads to large 
changes in the chemical composition of the disk, particu- 
larly in the cold and dense midplane layer. The effect of 
the accretion flow on a particular species depends on the 
accretion timcscale and the formation and destruction 
timescalcs of that species. For example, inside 10 AU, 
the gas-phase water abundances are not overly affected 
by the accretion flow, while ammonia is greatly enhanced 
and methanol, greatly depleted. 

Since the accretion flow affects mainly the midplane 
layers, the infrared molecular line emission is basically 
unchanged. Contrarily, turbulent mixing alters the 
chemical composition of the warm molecular layer, where 
the molecular line emission is generated. We thus con- 
trasted two models with and without turbulent mixing 
to observations of protoplanetary disks around classical 
T Tauri stars. We compared the modeled disk spectra 
and column densities of the upper layers of the disk to 
the observed spectra and the derived column densities. 
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In summary, a model with turbulent mixing repro- 
duces the observations of AA Tauri, DR Tauri and AS 205 
better than a model without mixing, although the dis- 
cussion of the column densities in the upper layers and 
of the spectra suggests a more efficient vertical mass 
transport mechanism. In the models considered here, 
turbulent mixing provides an efficient way of transport- 
ing material from the midplane into the transition layer. 
Due to the constant mixing velocity in each vertical col- 
umn, the turbulent timescales become longer than the 
chemical timescales at larger distances from the mid- 
plane. B ut these are the layers wh ere the disk winds 



tive (the transition layer coincides roughly with 1.5 disk 
scale heights) and drive large-scale channel flows. Hence, 
a combination of turbulence in the cold midplane layers 
and disk winds in the transition layer, both generated by 
the magneto-rotational instability, might account for the 
observations better than turbulent mixing or disk winds 
alone. In the latter case, the steady state assumption 
would have to be abandoned to allow for stronger disk 
winds. 

An interesting aspect for future studies will be the com- 
bination of radial and vertical transport processes. In 
this work, we showed that the radial accretion flow al- 
ters predominantly the midplane abundances, while tur- 
bulent mixing affects the midplane and warm molecular 
layers. We therefore expect that a combination of ra- 
dial and vertical motion will have a strong impact on the 
chemical abundances in the molecular layer, which will 
be reflected in the emission line spectrum. 
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Fig. 11. — Schematic view of (a) accretion, (b) mixing and (c) wind modeling in our model. 

APPENDIX 

A. PHYSICOCHEMICAL INTERACTIONS - NUMERICS 

Modeling the physical mass transport phenomena in combination with calculating the chemical evolution is the key 
aspect of this study. Here, we discuss in detail the numerical techniques used to compute the chemistry accounting 
for transport mechanisms in relation to the grid of our protoplanetary disk model. 

The general expression for the rate of change in the abundance of gas-phase species k in grid cell at radius s l 
and vertical distance z' J from the midplane is given by 

M _ ^ n k _ pi,3 _ r i,j _ - , • i,j ( A i \ 

k ~ At ~ k k k[cC fc,acc/mix/wind ' 

where Pk and Lk denote the production and loss terms due to gas-phase chemical reactions, hkice the adsorp- 
tion/desorption rate (> for net freeze-out, < for net evaporation) and fik , acc/mix/wind the rates arising from 
mass transport due to accretion, mixing or disk winds. In the case of molecular hydrogen, an additional (positive) 
term nn 2 is added to the right-hand side to account for the H2 formation on grains. Accordingly, in the case of atomic 
hydrogen, an additional term — 2tt,h 2 is required. 



A.l Radial transport by viscous accretion 

We assume that inward accretion occurs along constant streamlines z/h. (see Figure fTlTa) for an illustration of 
the coupling between the grid cells at different radii). The net accretion rate along constant streamlines is given in 
Equation (fT5|) . where the differentiation with respect to I needs to be translated into a differentiation with respect to 
s. Accordingly, on the discrete grid of the disk model, the net accretion rate for cell can be expressed as 

= ±r + e -p-^ 3 \ ( A2 ) 

where p and h represent the average mass and number density of the outer grid cells, weighed by their relative mass 
transfer into cell (i,j). The accretion velocity vl' %+1 is evaluated at the interface between the grid cells i and i + 1 

and is introduced to improve the accuracy of the finite calculation. The additional factor >•? j p 1 ^ arises from the 

i '-\-\ 1 2 

requirement of mass conservation. In the context of a steady disk model, all variables except and n k ' 3 ' 3 are 
constant with time. 



A. 2 Vertical transport by turbulent mixing 

Figure lllf b) demonstrates the grid-cell mass transfer by turbulent mixing in the diffusion-type model. Starting 
from Equation (fT8|) . we replace the differentials by finite differences over the vertical extent Az of the grid cells. 
Additionally, we introduce mean densities p and velocities v at the interface of the grid cells as before. The exchange 
rates of species k between grid cell (i, j) and the adjacent grid cells (i,j — 1) and + 1) then become 



' fc ,mix--^7=r — ; ^ ' ' \— "' ; "< : J • (A3) 

A. 3 Vertical transport by disk winds 

The model of the upward disk wind is shown in Figure [TTT c). In the finite case, the differential equation for the 
wind rate (Equation ([20]) ) of grid cell translates into 

. . v i,j-i,j . . , 7 ,».i,i+i . . 

U k, wind - Az i,j U k Az i,j n k ■ l A4 ^ 



In both the wind and the mixing scenarios, the only variable changing with time is nk- 
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